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1. Weight singularities and recurrence coefficients. 

The orthonormal polynomials p n (x) = 7 n x n + • • • related to the weight w satisfy the 
three-terms recurrence relation 

a n +iPn+i(x) = (x - b n )p n (x) - a n p n -i(x), (1) 

with aoP-i(x) = 0. 

Let w(x) > hold almost everywhere on the support [—1, 1], then one knows (since 1977) 
that the recurrence coefficients have limits a n — > 1/2 and b n — > when n — > oo (see for 
instance the survey in [Nev2]). 

Features of w can somehow be "read" in the sequence of the recurrence coefficients. This 
goes back to Stieltjes and is current practice in solid-state physics [LaG]. See [Ap] for the 
case of endpoints singularities. J. P. Gaspard once showed me a paper by C.Hodges [Ho] 
describing the influence of a mild interior singularity of the form 

w(x) ~ w(x ) + A(x - x ) 7 , x^x ,x>x , ^ 

~ w(xo) + B(xq — x) 7 , X^Xo,X<Xo, 

with x Q G (—1, 1),0 < w(x) < oo on (—1, 1) and 7 > (Van Hove singularity [Ho, Mart]) 
as 

a n — 1/2 ~ £n _7_1 cos(2?i(9o — if), b n ~ 2^?i _7_1 cos((2n + 1)6q —77), n — > 00 
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with xq = cos 6*o- I gave a lengthy proof of this in [Magi], showing how £ and r\ are related 
through Toeplitz determinants to the Szego function of w(cosO) sin# (the analytic function 
D(z) in \z\ < 1, with D(0) > 0, without zero in \z\ < 1, and such that the boundary values 
satisfy \D(e l0 )\ 2 = w(cos9) sin 9) by 

(^) 7+1 rh + D 

£ e «»7 _ V Z / e i(e -2arg(D(e ie 0)) _ B( ,-i-wy/2 

2tvw(xo) i 
Hodges argument is based on the continued fraction expansion of 

' x w(t) dt no 



z -t a? 

Z-Oq- 



when z = x + is for small e > and — 1 < x < 1, but the fastest explanation is probably 
related to inverse scattering techniques ( [VA] p. 117, [NV] and references therein): let us 
consider the function q n defined by the integral q n (z) = f_ 1 (z — t)~ 1 p n (t)w(t)dt when 
z 4- [-!>!]• Remark that f(z) = ^/Jj^q (z). From (1), a n+1 q n+1 (z) = (z - b n )q n (z) - 
a n q n -i(z) — y/JIoS n fi, (with cioq-i(z) = 0), which can be written as 

q n+1 (z)-2zq n (z)+q n - 1 (z) = q n+1 (z)-p(z)~ 1 q n (z)~p(z)[q n (z)-p(z)~ 1 q n - 1 (z)] =e n (z), n 



with e n {z) = (1 - 2a n+1 )q n+1 (z) - 2b n q n (z) + (1 - 2a n )q n - 1 (z) - 2 y /p ~5 rLj0 , and where 
p(z) is the determination of (z + \J z 2 — l)/2 such that \p(z)\ > 1 when z £ [—1, 1]. After a 
simple summation, p(z)~ N (q N+1 (z) - p(z)~ 1 q N (z)) - p(z)q (z) = Y^o p(z)~ n e n (z). While 
z £ [—1, 1] and when N — > oo, p(z)~ N and gjv(^) — > 0, so 

/ x = fjz) _ £n(z) _ 2 y / ^o v 2 ^ (2a n - l)(g n _i(^) + /o(z)g n (z)) + 2b n q n (z) 
Q0[Z) ~ VJ^o " ~ V ^^) n+1 ~ o ^ 



n+1 



From the known asymptotic behaviour ~ , 7 . - — , ;t~tt^7 7 f° r large n [Ba], we 

(p(z) - p(z) 1 )p(z) n 

have a first approximate expansion of {2tt)~ 1 ' 2 D{p{z)~ 1 ){p{z) — p(z)~ 1 )qo(z) as a series 
of negative powers of p(z): 

Under sufficiently strong conditions , this remains valid when z = x+ie, e — > 0, £ > 0, — 1 < 
a; = cos(#) < 1, p(z) — > e (with < # < 7r), so that 2a n — 1 and b n behave like the 
Fourier coefficients of order 2n and 2n + 1 of a function whose singularities on < 6> < 7r 
are related to the singularities of iu(cos6>) sin6>. Singularities of type \6 — 6>o| 7 correspond 
indeed to n -1-7 cos(n6>o+ const.) behaviour in the n th Fourier coefficient. 
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This analysis is no more valid for stronger singularity w(xq) = or oo because the 
approximations done near xq are no more valid, subtle important effects take place in 
neighbourhoods of length about 1/n of the singular points, see the famous [NevIII]. 

It is therefore not useless to have a close look at the simplest orthogonal polynomials 
related to weights with interior singularities. Special singular positions have been worked 
(in sieved polynomials theory etc.), but here is something related to an arbitrary 
position: 

2. Freud's equations for the simplest generalized Jacobi polynomials. 

Let 



W(x) = 5(1 - X) a (x - XY(\ +xf XE [-1, Xq] , 

= A(l - x) a (x - x ) 7 (l + xf xe [x Q , 1] , 



e ( 3 ) 



with — 1 < xq < 1, A and B > 0, a, (3 and 7 > — 1. 

The case A = £?, a = (3, xq = has the simple solution deduced from Jacobi polynomials 
a 2 n = (n + 2a + 7odd(n))(rc + 7odd(n))/[(2w + 2a + 7 + l)(2n + 2a + 7 - l)],b n = 0, 
where odd (n) = (1 — (— l) n )/2. When n — > 00, this case shows the asymptotic behaviour 

2 An 

The 0(l/n) term is definitely related to the |x| 7 behaviour of the weight near 0, as shown 
by P. Nevai (Theorem 4 of Section 7 of [Nev]): if w is even on [—1,1], with w(x)\x\~ 1 
positive and continuously derivable on (—1,1), then a n = 1/2 — 7(— l) n /(4n) + o(l/n) 
when n — > 00. 

We now try to investigate the recurrence coefficients a n , b n when the weight is (3). This 
weight is a semi- classical weight, as w'/w is the same rational function almost everywhere 
on the support [—1,1] of w. 

Semi-classical orthogonal polynomials have a rich differential structure, according to a 
theory going as far as Laguerre [BeR, GaN, Lag, Mag2, Mag3, Sho]. 

Freud [Fr] showed how to deduce recurrence coefficients asymptotics from special 
identities. For a general semi-classical weight satisfying w'(x)/w(x) = 2V(x)/W(x) with 
W{x)w{x) — > when x tends to any endpoint of the support S of w (Shohat's conditions 
[Sho]), we find these identities (Freud's equations) by expanding 



0= [W(x)w(x)p n {x)p n - k (x)]' dx 
Js 

= / W(x)w(x)p' n (x)p n -k(x) dx + / W(x)w(x)p n (x)p' n _ k (x)dx+ 
Js Js 

+ / W' (x)w(x)p n (x)p n -k{x) dx + / 2V(x)w(x)p n (x)p n - k (x) dx 
Js Js 

for k = 0,1, remarking that any integral f s P(x)w(x)p n (x)p n -k(x)dx where P is a 
polynomial, is an expression involving a n , b n , a n ±i, b n ±i, etc. according to k and the 

degree of P, and using p' n = np n -i/a n + (b H h6 n _i - nb n -i)p n - 2 / (a n -ia n ) H (see 

[BeR]). 
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Reduction to even measure: one considers the orthonormal polynomials {p n } with respect 
to the even weight 

w(x) = 2\x\w(2x 2 - 1), for - 1 < x < 1. 

Then, p2 n (x) = Pn(2x 2 — 1), and one recovers the recurrence relation for the p^s by 
contracting the recurrence relation for the p n 's: 

d n+1 p n+1 (x) = xp n (x) - a n p n -i(x) =>• 

a2n+ia2n+2P2n+2(x) = (x 2 - d\ n - 0?, n+1 )p2 n (x) - a 2 n-ia 2 nP2n-2 (x) , 

so: 

2a 2n -ia 2n , b n = -1 + 2d\ n + 2dl n+1 (A) 



a r 



This allows to work with the single sequence {a n } instead of the two sequences {a n }, {b n }. 
Here, 

w(x) = 2\x\w(2x 2 - 1) = B\x\ 2f3+1 (xl - x 2 y(l - x 2 ) a for \x\ < \x \, 

= A\x\ 2f3+1 (x 2 - 5-^(1 - x 2 ) a for |x | < \x\ < 1, 

where x is the positive root of 2x% - 1 = x , A = 2 a + f3 +' y+1 A, B = 2 a + f3 +' y+1 B. 

So, W(x) = x(x 2 - xl)(x 2 - 1) = x 5 - (£q + l)x 3 + XqX and 
2V(x) = (2a + 2(3 + 2^ + l)x 4 - [2ax 2 + (2(3 + l)(x 2 + 1) + 2 7 ]x 2 + (2(3 + l)x 2 . 

The equations for the a n 's now follow from Freud's method for even weights, expanding 
Ww' = 2Vw as 

f W(x) „, f W(x) f ( W(x)\ 

/ W(x)p' n (x)p n - 1 (x)+ W(x)p n (x)p' n _ 1 (x)+ W (x)p n (x)p n - 1 (x) + 

Js x j s x j s y X J 

f 2V(x) „, _ , _ , . n 
+ / w(x)p n (x)p n -i(x) = 0, 

using f s w(x)pi(x)pj(x)dx = dij, the recurrence relations (1) giving f s xw(x)p n (x)p n -i(x)dx 
a n , f s x 2 w(x)(p n (x)) 2 dx = d\ + d 2 l+1 , f s w(x)x~ 1 p n (x)p n -i(x)dx = odd(n)/a n , etc., and 

Pn = —Pn-l H = = = Pn-3 + 

a n a n -2d n -ia n 
ndl^d 2 ^ - 2(d 2 n _ 3 + d 2 n _ 2 + dl^) Ei" 1 ~ a l + 2 Ei _1 (at + 2 ~ a l~ a l-i) ~ 

H Z Z Z Z Z Pn-b H , 

fln-4«n-3 a n-2«n-lfln 

one finally finds 

2(n + a + (3 + 7 + 2)^(^_ 1 + d 2 n + d 2 n+1 ) - 2[ax 2 + (n + (3 + l)(x 2 + 1) + 7 ]a 2 + 

n— 1 n— 1 

+2(2a 2 - x 2 - 1) ^ a 2 + nx 2 - 2d 2 n d 2 n _ 1 + 2 ^(^ 4 + 2d 2 d 2 _ 1 ) + (2(3 + l)x 2 odd (n) = 0, 
i=i i=i 

n = l,2,... (5) 



4 



Simplest generalized Jacobi polynomials. 

(o = 0). 

We see how any a n can be computed from the value of Si, which is the degree of freedom 
reflecting that the same equations (5) hold for any choice of A and B in the weight w. 
Actually Si is linked to the ratio A/B by 



a i 



&2 

B L, - \x\W+\5% - x 2 )-y(l - x 2 ) a dx + A f , \x\ 2 ^{x 2 - - x 2 ) a dx 



Bf ]x]< - o \x\ 2 P+ 1 {xl-x 2 )~<{l-x 2 ) a dx + Aj\ x>io \x\ 2 ^{x 2 -x 2 p(l -x 2 ) a dx 

(6) 



Numerical experiments show that S 2 , as,... can be computed in a stable way from Si 
simply by considering (5) as an equation for S n +i when Si, . . . , d n are known: 

~ 2 _ ax 2 + (n + /3 + l)(x 2 + l)+<y . Ei" 1 ^ , 
a -+i" N N 

, 2(xg + l) Er 1 4-^-2Ei" 1 («t + 2SlSl_ 1 )-(2/3 + l)x 2 odd(n) _ a 2 _i . 2 . 2 



2iVS n AT 



a n a n-l 



(7) 



with A r = n + a + /3 + 7 + 2, n = l,2, .... 
3. Asymptotic estimates. 

Some people call these tricks "special refinements" ; 
others call them "kludges" . 
D.E. Knuth 

Putting an almost constant a 2 ~ S 2 in (5) gives two possible asymptotic matches a 2 = 1/4 
and S 2 — Xq/4, corresponding to weights with support [—1,1] and [— xo, xo] (the latter when 
A = 0). More complicated behaviours are expected to hold when x is complex [GaN] and 
one should be able to establish correct asymptotic behaviours from (5) alone, but this 
has not yet been achieved (the answer to [GaN] was given in [N] with other techniques). 
Even when one knows that a 2 — > 1/4, some amount of guesswork will still be needed. Let 

a 2 = - + y n . We know from Szegd's theory ([Sz] chap. 12) that Ei a \ = ^ + £ + z n 
and Ei~ (®k + 2S|S|_ X ) = — + rj + u n with z n and u n — > (See also [Nev2] p. 91). So, 

n 

1 axl + (n + f3 + + 1) +7 n 4 + ^ + ^ n 
2/n+i = -7 H 77 2 



4 A" A" 

3n 

16 



2(xg + l)(- + e + ^)-nxg-2(— + 7/ + u n )-(2^+l)xgodd(n) 7+2/n-i 1 



-^-Vn-Vn-1- 
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We compute now £ and r\: from the Szego theory, let 4> n {z) = k n z n +k' n z n ~ 2 +k" n z n ~ 4: -\ 

be the orthonormal polynomials on the unit circle with respect to | sm9\w(cos9) = 
C(6)\cos6\ 2 ( 3+1 \cos 2 6-cos 2 (6o/2)\'y\sm6\ 2a+ \ with C{9) = A on \9\ < 9 /2 and 
\9 — 7r j < 9q/2 and C(9) = B elsewhere. The Szego function D(z) whose boundary 
values must be \D(e l9 )\ = ^w(cos 9)\ sin 6*| is found by inspection to be 

D(z) = R-\l - z 2 r +1 / 2 (l + zy+^il - e- i0 °z 2 y' 2+iX (l - e^z 2 yl 2 - lX , 

with k = 2 a +P+~<+ 1 B( eo - 7T ')/( 2 ^A- 9o /( 2 ^ and A = (In)' 1 \og(B/A). We find the limit 
values of k' n and k" n from the expansion of 1/D: 

fc'n/ a — P + l x o + 2A sin 9q, 
Kn/K n -> (a-/3) 2 /2 + (a + /3 + l)/2 + (a-/3)(7x + 2Asin^ ) + 

+(7(7 + 2)/4 - A 2 ) cos(26> ) + A( 7 + 1) sin(26> ) + 7 2 /4 + A 2 , 



used in 

Pn(x) = 



Pnjx) 
In 



= X 



'n-\ 



x n ~ 2 + 



'n-\ 



'n-\ 



,fc=l 



vfc=l 



2 n (k 2n + 4> 2n (0) 
T n (x) k' T n _ 2 (x) k" T n _ 4 (x) 



2 n-l 



H 2 n ~ 1 

x 1 



K 



2 n-l 



'n-2 



- E4 -2 £ 



~2 ~2 
a k a k+l 



vfc=l 



_ n - 2 n(n - 3)/2 - (n - 2)^//t + k/^/k 4 

r *- / Jb Jb \~ t Jb 



16 



x n - 4 /2 + 



(z + z' 1 
So, 



2x). 



n-l 



£ = lim \^ a 2 k - n/4 = -k'/(4k) = -(a - f3 + <yx + 2A sin O ) /4, 



(8) 



7/ = lim^-.oo Ei ^Ofc + 2afa 2 _ 1 ) - 3n/16 = {{k' / k) 2 - 4k' /k - 2k"/k)/16. 
The equation for the y n 's reduces to 



y n+ i - 2x y n + y n -i = 



(x + !)(/? + l/2)(-l) n + 2(x + 2)z w + y n _! - 4n n 

AT 



16y n (2x + l)/4 - A sin 9 + (xg + 1) Q3 + 1/2) (-1)" + 2(x + 3)z n - 4u n 
l + 4y n 4iV 

4(2xo + l)vl 
l + 4y n 



(9) 
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where z n = - J2n Vk and u n = -J/n-i/2 - X^( 3 W 2 + vt + ^VkVk-l)- 

The form (9) should give hints on the behaviour of y n when n — > oo. No proof will 
be attempted here, only reasonable asymptotic matching and numerical checks. Use of 
Painleve-like differential equations in 9q is another method of investigation which could be 
used in the future (see [Mag3]). 

The right-hand side is small, even with respect to the j/'s, so y n +i — 2xoy n + Un-i is 
small, and this suggests a exp(±m#o) behaviour somewhere. However, I still don't have a 
tight proof that y n , z n and u n are 0(l/n). Assuming y n = Ki(—l) n /n + K 2 e /n^ 2 + 
X 3 e~ /n^ 3 , matching the two sides gives K\ = —(3/2 — 1/4, (2,3 = 1 ± 2i\. Numerical 
checks have been performed on the form 

a 2 - 1/4 = y n = -(J3 + l/2)(-l)7(2ra) + Kcos(n6> - 2Alogn - ip)/n + o(l/n) (10) 

where K and ip are unknown functions of a, /?, 7, A and xq. Given a, 7 and 
xo, the algorithm first performs (7) with several trial starting values a\ and computes 
the corresponding A from (8), allowing the determination of the coefficients in (6). It 
is then possible to run (7) for a requested value of A, and to estimate K and ip in 
(10) from numerical values of for large n (up to the 10000-100000 range). Very 
satisfactory empirical formulas for K and ip appear to be K = (7 2 /4 + A 2 ) 1 / 2 sin(6*o/2) and 
<p>= (a + l+ 7 /2)7r-(« + /3 + 7 + l/2)6'o + 2Alog(2sin6'o)-2argr(7/2 + zA)-arg(7/2 + zA). 
Whence, from (4): a n - 1/2 ~ y 2n -i + V2n, K ~ 2(y 2 n + 2ten+i), the 

Conjecture. The recurrence coefficients related to the simplest generalized Jacobi weight 
(3) satisfy 

1 M 

a n = cos [2n$o — 2A log(4n sin ^ ) — + o 

b n = - — cos \(2n + l)6»o - 2A log(4n sin 6 )-$] + o 

n 

when n -> 00, where x = cos9 , < 9 < n, A = \og(B /A)/(2ir), M = ^(7 2 / 4 + 
A 2 ) 1 /2 sin ^ 0) $ = ( Ck + 7 /2) 7r - ( a + /3 + 7 )6» -2argr( 7 /2 + zA) -arg( 7 /2 + zA). 

Here is a sample of the numerical check: K and </? are extracted from the form (10) on a 
sample of a 2 , with n going up to 500000. The values of K and A (from (8)) are quite stable, 
but things are not so easy with ip (the "phase" column). Finally, the line "check" contains 
the computed values of K and <p, and a new check of A through Turan determinant weight 
function reconstruction yielding A and B. 

Script started on Wed Mar 2 17:38:49 1994 
/ul8/grpanma/magnus/jacobi@uxl2 [1] #ls -1 
total 1224 

-rwxr-xr-x 1 magnus 473298 Mar 2 17:30 a. out 
-rw-r— r— 1 magnus 2226 Feb 16 15:24 alnrel.f 

-rwxr-xr-x 1 magnus 18 Jan 3 14:15 batchsub 

-rw-r 1 magnus 1183 Feb 16 08:58 c81gmc.fz 
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-rw-r— 





1 


magnus 


2087 


Feb 


16 


08: 


58 


c91gmc . f 


-rw-r- 


-r — 


1 


magnus 


259 


Feb 


16 


15: 


25 


carg.f 


-rw-r- 





1 


magnus 


21751 


Feb 


16 


16: 


41 


clngam.f 


-rw-r- 


-r — 


1 


magnus 


14412 


Feb 


16 


16: 


42 


clngam. o 


-rw-r- 





1 


magnus 


2645 


Feb 


16 


15: 


39 


clngamold.f 


-rw-r- 





1 


magnus 


975 


Feb 


16 


08: 


58 


clnrel . f 


-rw-r- 


-r — 


1 


magnus 


141 


Feb 


16 


16: 


35 


cmplxt . f 


-rw-r- 


-r — 


1 


magnus 


1061 


Feb 


16 


15: 


52 


csevl . f 


-rw-r- 





1 


magnus 


471 


Feb 


16 


08: 


58 


entsrc . f z 


-rw-r- 


-r — 


1 


magnus 


7679 


Mar 


2 


17: 


29 


gj acobi.f 


-rw-r- 


-r — 


1 


magnus 


14554 


Mar 


2 


17: 


30 


gj acobi . o 


-rw-r- 


-r — 


1 


magnus 


6394 


Feb 


9 


09: 


07 


gj acobi4.f 


-rw-r- 





1 


magnus 


4153 


Feb 


16 


08: 


59 


indexfn.z 


-rw-r- 


-r — 


1 


magnus 


861 


Feb 


16 


15: 


52 


inits . f 


-rw-r- 


-r — 


1 


magnus 


167 


Feb 


16 


16: 


46 


lgam.f 


-rw-r- 


-r — 


1 


magnus 


1325 


Feb 


16 


16: 


46 


lgam. o 


-rwxr- 


xr-x 


1 


magnus 


452258 


Mar 


2 


13: 


50 


lgam. out 


-rw-r- 




1 


magnus 


10873 


Feb 


16 


08: 


58 


macharb . f 


-rw-r- 




1 


magnus 


4759 


Feb 


16 


08: 


58 


rlmacha. f z 


-rw-r- 




1 


magnus 


8318 


Feb 


16 


08: 


58 


rlmachb . f 


-rw-r- 




1 


magnus 


2800 


Feb 


16 


08: 


58 


seterr . f 


-rw-r- 




1 


magnus 


319 


Feb 


16 


08: 


58 


seteru.f 


/ul8/grpanma/magnus/ j acobi@uxl2 


[2] 


#cat gjacobi.f 



gj acobi. for 

program gj acobi 
simplest generalized Jacobi 



w(x) = B (1-x) "alpha 
= A (1-x) "alpha 

-> w tilde (y)= B tilde 

-> w tilde (y)= A tilde 



(1+x) "beta 
(1+x) "beta 



(xO-x) "gamma 
(x-xO) "gamma 



|y| "(2beta+l) 
|y| "(2beta+l) 



(y0~y) "gamma 
(y~y0) "gamma 



-Kx<x0 
x0<x<l 

(l-y~2) "alpha 

-y0< |y | <y0 
(l-y~2) "alpha 

y0<|y|<l 



x=2y"2-l 

the program produces squares of recurrence coefficients atilde n 

of atilde n+1 ptilde n+1 

where the ptilde s are the 
input 



(y) = y ptilde n (y) - atilde n ptilde n-1 (y) 
orthonormal polynomials related to w tilde 



first, period of display of atilde n, number of values of y 

where to compute polynomials, and interval (this latter part 
is alpha test, enter 0) 

then: alpha, beta, gamma, yO , n max 

then, after preliminary calculations have been done, 
lambda = log(B/A) /(2 pi) 

integer oddn,n,nmax, i , signg,recomp, period, npol 
double precision al,be,ga,yO,atn,atnpl,pi,aux,shalf ,cth,sth, 
& at2 , at2pl , at2ml , one , half , 

& thO, x, fn, cnt,snt,sl,s2,snpl,cnpl,xx,xy,yy,laml,lam2 

& ,auxl,at,bt,ba(4) ,at21(4) ,cl,c2,c3,thx0,lam0,phi0,phi 

& ,phil,lsth 
or, if you can afford it (useful if | lambda |>5) : 
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c real*16 al ,be ,ga,yO , atn, atnpl ,pi , aux, shalf , cth, sth, 
c & at2 , at2pl , at2ml , one , half , 

c & thO, x, fn, cnt,snt,sl,s2,snpl,cnpl,xx,xy,yy,laml,lam2 

c & , auxl , at ,bt ,ba(4) ,at21(4) ,cl,c2,c3,thx0,lam0,phi0,phi 

c & ,phil,lsth 

double precision xv (502) ,polv (502,3) ,dtur 
c complex gamma function from W. Fullerton's fnlib 

c available by anonymous ftp on netlib.att.com cd/netlib/fn 

complex clngam,z 
one=l 

half =one/2 
pi=4*atan(one) 

print *,' period of display of coef f . , number of abs . , extr. abs.? ' 

read * .period, npol ,xvl ,xv2 
c recurrence coefficients will be printed with step period, 
c polynomials and weight will be computed at xvl,...,xv2 

npol=min(npol , 500) 
1 print *,' input alpha, beta, gamma, xtilde0,nmax (stop if nmax<=0) ' 

read *,al,be,ga,y0,nmax 

if (nmax . le . 0) stop 

print ' (4f 14.8, i8,f 9.4) ' ,al,be,ga,y0,nmax 
signg=l 

if (ga. It .0.0) signg=-l 
shalf=sqrt(l-y0*y0) 
cth=2*y0*y0-l 
sth=2*shalf *y0 
th0=acos (cth) 

phi0=pi*(2*al+2+ga)/2 - (al+be+ga+half ) *th0 
lsth=log(2*sth) 
xv(l)=y0/2 
xv(2)=(l+y0)/2 
do 100 i=l,npol 
100 xv(i+2)=xvl+(i-l)*(xv2-xvl)/(npol-0.999) 
itba=l 
at2=y0*y0 
at21(l)=at2 

c 3 satisfactory trial values of atilde square(l) are needed to 
c establish the homographic relation atilde square(l)=cl+c2B/(A+c3B) 
11 if (itba.eq. Dgoto 12 
if (itba.le.3)then 
if (at21 (itba) . le . 0) at21 (itba) =0 . 0001 
if (at21(itba) .ge. I)at21 (itba) =0.9999 
at21 (itba)=0 . 99*at21 (itba) +0 . 01*y0*y0 
c to be tried until all the atilde square (n) are positive 
at2=at21(itba) 
else 

print *,' ? lambda ( stop if > 999)' 

read *, lamO 

if (lamO.gt. 999) goto 1 

ba(4)=exp(2*pi*lam0) 

at2=cl+c2/(c3+l/ba(4)) 

at21(4)=at2 
endif 
12 cnt=cth 
snt=sth 
sl=0 
s2=0 
at2ml=0 
oddn=l 
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if (itba.ge.4) 

feprint *,' n atilde2(n) K phase lambda lambda' 

if (itba.ge.4) print ' (i8 ,f 14 . 10) ' , 1 , at 2 
atn=sqrt (at2) 
do 13 i=l,npol+2 
polv(i , 1)=1 
13 polv(i,2)=xv(i)/atn 

main engine for recurrence coefficients 
do 10 n=l,nmax 

aux=al*y0*y0+(n+be+l)*(y0*y0+l)+ga+( (y0*y0+l) /at2 -2 )*sl 
& -(n*y0*y0)/(2*at2)+at2ml-s2/at2 - (be+half ) *y0*y0*oddn/at2 
at2pl=aux/ (n+al+be+ga+2) -at2ml-at2 
updating the two sums 
sl=sl+at2 

s2=s2+at2*( at2+2*at2ml ) 
snpl=snt*cth+cnt*sth 
cnpl=cnt*cth-snt*sth 
check asymptotic formula when needed to display 

if ((n. It. 2) .or. (at2pl.lt.O) .or. (mod(n+l .period) .eq.O) .or. 
& (n.eq.nmax)) then 

aux=(n+l) * (at2pl-half /2) +(oddn-half ) * (be+half) 
aux should be about K cos( (n+l)theta0 -phi ) 
with phi slowly varying with n 

auxl=n* (at2-half /2) - (oddn-half ) * (be+half) 
K cos (phi) : 

xx=(auxl*snpl-aux*snt) / sth 
K sin(phi) : 

y y = ( aux * cnt - aux 1 * cnp 1 ) / s t h 
phi (between -pi and pi ) 
phi=atan2 (yy , xx) 

K 

xx=sqrt (xx*xx+yy*yy) 
2 ways to estimate lambda 

laml=(-4*sl+n+l-al+be-ga*cth)/(2*sth) 

lam2=half * (-16*s2+3* (n+1) -5*al+3*be-l-ga* (2*cth*cth+4*cth-l) ) / 
& (2*sth*(cth+2)) 
fn=n+l 
phi -2 lambda log n 

f n=phi-2*laml*log(f n) 
fn=fn/(2*pi) 
if n=f n 
f n=f n-if n 

if (f n . It . -half ) fn=f n+1 
if(fn.gt. half)fn=fn-l 
phi=2*pi*f n 
Theta n(x tilde 0) : 

thx0=4*sl+2* (n+al+be+ga) *at2+2* (n+2+al+be+ga) *at2pl 
& -2*(n+be+ga+l-ga*y0*y0) 
thx0=thx0*y0 
if (itba.ge.4) 

feprint ' (i8,f 14.10, 4f 10.5) ' ,n+l,at2pl,xx,phi,laml,lam2 
endif 

if (at2pl.le.0)then 
if (itba.ge.4) 

& print *,' atildesquare <0, try a less extreme lambda' 

goto 11 

endif 

values of orthogonal polynomial of degree n+1 
atnpl=sqrt (at2pl) 
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do 25 i=l,npol+2 
25 polv(i,3)=( xv(i)*polv(i,2)-atn*polv(i,l) )/atnpl 
prepare next step 

if (n . It . nmax) then 
snt=snpl 
cnt=cnpl 
oddn=l-oddn 
at2ml=at2 
at2=at2pl 
atn=atnpl 
do 251 i=l,npol+2 
polv(i , l)=polv(i ,2) 
251 polv(i,2)=polv(i,3) 

endif 
10 continue 

if (itba. ge . 4) then 
check : weight function reconstruction 
w(x)= 2 sqrt(l-x~2) /(pi * dtur) 
do 20 i=l,npol+2 
x=xv(i) 

dtur=polv (i , 2) **2-atnpl*polv (i , 1) *polv (i , 3) /atn 
fn=l/(dtur*pi* x**(2*be+l) * (abs (y0*y0-x*x) ) **ga 
& *(l-x*x)**(al-half )) 

if (i.eq. l)bt=fn 
if (i . eq. 2) at=f n 

if (i.gt.2)print ' (lp,5el5.7) ' ,x,fn, (polv(i, j) , j=l,3) 
20 continue 

print '(" Atilde=" ,lp,el2.4, " Btilde=" ,el2.4, " B/A=", 
& el2.4) ' ,at,bt,bt/at 
phil=phi0+2*laml*lsth 
z=cmplx(ga/2 , laml) 
z=clngam(z) 

phil=phil-2*aimag(z)-atan(2*laml/ga) 

phil=phil+(l-signg) *pi/2 

phil=phil/(2*pi) 
iphil=phil 
phil=phil-iphil 
if (phi 1 . gt . half ) phi l=phi 1- 1 
if (phil . It . -half ) phi l=phi 1+1 
phil=phil*2*pi 

print '(" check: " ,15x,3f 10 . 5 , lp, 2el0 . 3) ' , 
& shalf *sqrt(ga**2+4*laml**2)*half , 

& phil, half *log(bt/at)/pi, at, bt 

endif 

estimate of B/A = exp(2pi lambda) 

ba(min(itba,4))=exp(2*pi*laml) 

if (itba.lt.4)print ' (i3,f 15.8) ' , itba, at21 (itba) 

itba=itba+l 

recomp=0 

if (itba. eq. 4)recomp=l 
recompute the homographic relation if new extreme data have been found 
if (itba. gt. 4) then 

if (at21 (4) . It . at21 (2) ) then 
recomp=l 
at21(2)=at21(4) 
ba(2)=ba(4) 
endif 

if (at21(4) .gt.at21(3))then 
recomp=l 
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at21(3)=at21(4) 
ba(3)=ba(4) 
endif 
endif 

if (recomp . eq. 1) then 
xx=(at21(l)-at21(2))/(l/ba(l)-l/ba(2)) 
yy=(at21(l)-at21(3))/(l/ba(l)-l/ba(3)) 
c3=(l/ba(3)-l/ba(2))/(xx/yy-l)-l/ba(2) 

c2=(at21(l)-at21(2))*(l/ba(l)+c3)*(l/ba(2)+c3)/(l/ba(2)-l/ba(l)) 
cl=at21(l)-c2/(l/ba(l)+c3) 

print ' (" (" ,f 15.6, " A+",fl5.6," B)/(A+ ",fl5.6," B)")', 
& cl,cl*c3+c2,c3 

print '('' al tilde square min,max= ' ' ,2f 10.5) ' ,cl,cl+c2/c3 
endif 

if (itba.eq.2)at21(2)=0 
if (itba.eq.3)at21(3)=l 
goto 11 
end 

/ul8/ grpanma/magnus/ j acobi@uxl2 
period of display of coef f . , 
50000 

input alpha, beta, gamma, xtilde0,nmax 
0.1 -0.66 -0.9 0.8 500000 

0.10000000 -0.66000000 -0.90000000 
1 0.64000000 2 0.48642906 3 

( 0.666359 A+ 0.682666 B) / (A+ 

al tilde square min,max= 0.66636 0.48607 
? lambda ( stop if > 999) 1 



[3] #a.out 
number of abs . 



extr . abs . ? 



(stop if nmax<=0) 



0.80000000 500000 
0.66612028 
1.404457 B) 



n 


atilde2(n) 




K 


phase 


lambda lambda 


1 


0. 


.4863105269 


















2 


0. 


.1072928956 





.56893 


0. 


63299 


-0, 


.23606 


-0, 


.11591 


50000 


0. 


.2500104918 





.65794 


2. 


04846 


1. 


.00002 


1, 


.00003 


100000 


0. 


.2499993374 





.65795 


2. 


04921 


0. 


.99999 


1, 


.00000 


150000 


0. 


.2500017530 





.65794 


2. 


04881 


1. 


.00001 


1, 


.00001 


200000 


0. 


.2499986110 





.65795 


2. 


04911 


1, 


.00000 


1. 


.00000 


250000 


0. 


.2500025306 





.65795 


2. 


04887 


1. 


.00001 


1. 


.00001 


300000 


0. 


.2499980788 





.65795 


2. 


04903 


1, 


.00000 


1. 


,00000 


350000 


0. 


.2500018867 





.65796 


2. 


04897 


1, 


.00000 


1. 


.00001 


400000 


0. 


.2499994423 





.65795 


2. 


04895 


1. 


.00000 


1. 


,00001 


450000 


0. 


,2499999531 





.65795 


2. 


04904 


1, 


.00000 


1. 


,00000 


500000 


0. 


.2500011150 





.65795 


2. 


04892 


1, 


.00001 


1. 


.00001 


500001 


0. 


.2500009765 





.65795 


2. 


04898 


1, 


.00000 


1 


.00001 


check: 









.65795 


2. 


04898 


1. 


.00000 


6.315E-05 3.382E-02 


0. 


,666359 A+ 


0.682666 


B)/(A+ 




1.404457 B) 


al tilde 


square min,max 




0.66636 


0.48607 








? lambda 


( stop if > 


999) 


2 












n 


atilde2(n) 




K 


phase 


lambda lambda 


1 


0. 


.6663578695 


















2 


0. 


.0057517038 





.88852 


1. 


15619 


-0 


.61116 


-0, 


.49518 


50000 


0. 


.2500181130 


1 


.22984 


-2. 


15364 


-1, 


.99944 


-1. 


.99940 


100000 


0. 


.2500124101 


1 


.22976 


-2. 


15333 


-1, 


.99945 


-1 


.99943 


150000 


0. 


.2500050485 


1 


.22973 


-2. 


15338 


-1, 


.99945 


-1. 


.99944 


200000 


0. 


.2499972423 


1 


.22969 


-2. 


15318 


-1, 


.99946 


-1. 


.99946 


250000 


0. 


.2499956982 


1 


.22970 


-2. 


15294 


-1, 


.99947 


-1. 


.99947 


300000 


0. 


.2500015654 


1 


.22970 


-2. 


15304 


-1, 


.99947 


-1. 


.99946 


350000 


0. 


.2500034533 


1 


.22971 


-2. 


15326 


-1, 


.99946 


-1 


.99945 


400000 


0. 


,2499985243 


1 


.22969 


-2. 


15320 


-1, 


.99946 


-1. 


,99946 


450000 


0. 


.2499982891 


1 


.22970 


-2. 


15304 


-1, 


.99947 


-1. 


,99946 


500000 


0. 


.2500023390 


1 


.22969 


-2. 


15318 


-1, 


.99946 


-1. 


.99946 
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500001 0.2499993553 1.22970 -2.15306 -1. 
check: 1.22969 -2.15322 -1. 

0.666359 A+ 0.682666 B)/(A+ 

al tilde square min,max= 0.66636 0.48607 
? lambda ( stop if > 999) 3 



99947 -1.99946 
99946 4.756E-02 
1.404457 B) 



1.664E-07 



n 


atilde2(n) 




K 


phase 


lambda lambda 


1 


0. 


.4860711314 


















2 


0. 


.1073891907 


0. 


56857 


0. 


.63221 


-0. 


23556 


-0. 


.11549 


50000 


0. 


,2499658427 


1. 


90240 


1. 


.71797 


3. 


.13855 


3. 


13861 


100000 


0. 


,2499817770 


1. 


90236 


1. 


.71765 


3. 


13858 


3. 


13861 


150000 


0. 


,2499992879 


1. 


90235 


1. 


.71679 


3. 


13862 


3. 


,13864 


200000 


0. 


,2500085823 


1. 


90243 


1. 


.71702 


3. 


13861 


3. 


13864 


250000 


0. 


.2499927104 


1. 


90240 


1. 


.71719 


3. 


13861 


3. 


.13862 


300000 


0. 


,2500063733 


1. 


90242 


1. 


.71677 


3. 


.13862 


3. 


13864 


350000 


0. 


,2499951650 


1. 


90239 


1. 


.71715 


3. 


13861 


3. 


13862 


400000 


0. 


.2500046810 


1. 


90242 


1. 


.71677 


3. 


13862 


3. 


.13864 


450000 


0. 


.2499960453 


1. 


90241 


1. 


.71706 


3. 


13861 


3. 


13862 


500000 


0. 


.2500039648 


1. 


90243 


1. 


.71682 


3. 


13862 


3. 


.13863 


500001 


0. 


,2500009249 


1. 


90243 


1 


.71704 


3. 


.13862 


3. 


13863 


check: 






1. 


90243 


1. 


.71683 


3. 


13864 


9.229E-11 3.386E-02 





.666359 A+ 


0.682666 


B)/(A+ 




1.404457 B) 



al 

? 



tilde square min,max= 0.66636 
lambda ( stop if > 999) -4 



0.48607 



n 


atilde2(n) 




K 


phase 


lambda lambda 


1 


0. 


.6663587554 


















2 


0. 


.0057510998 


0, 


.88853 


1. 


15619 


-0. 


.61116 


-0. 


.49518 


50000 


0. 


.2500178596 


2. 


.17836 


-0. 


90688 


-3. 


.60233 


-3. 


.60219 


100000 


0. 


,2499851215 


2, 


.17817 


-0. 


90842 


-3. 


.60225 


-3, 


.60221 


150000 


0. 


,2499982856 


2, 


.17817 


-0. 


90898 


-3. 


.60222 


-3. 


.60219 


200000 


0. 


,2499966056 


2, 


.17815 


-0. 


90894 


-3. 


.60222 


-3, 


,60220 


250000 


0. 


,2499928901 


2, 


.17814 


-0. 


90874 


-3. 


.60223 


-3, 


.60221 


300000 


0. 


,2499936079 


2, 


.17816 


-0. 


90855 


-3. 


.60224 


-3, 


.60222 


350000 


0. 


,2500002389 


2, 


.17818 


-0. 


90867 


-3. 


.60223 


-3, 


.60221 


400000 


0. 


,2500054650 


2, 


.17815 


-0. 


90898 


-3. 


.60222 


-3, 


.60220 


450000 


0. 


,2500020824 


2, 


.17816 


-0. 


90913 


-3. 


.60221 


-3, 


.60220 


500000 


0. 


.2499962590 


2, 


.17814 


-0. 


90891 


-3. 


.60222 


-3, 


.60221 


500001 


0. 


,2500006090 


2, 


,17815 


-0. 


90911 


-3. 


.60221 


-3, 


.60220 


check: 






2, 


.17813 


-0. 


90913 


-3. 


.60221 


4.756E-02 7.042E-12 





.666359 A+ 


( 


). 682666 B)/(A+ 




1.404457 B) 


al tilde 


square min,max= 


0.66636 


0.48607 








? lambda 


( stop if > 


999) 


5 












n 


atilde2(n) 




K 


phase 


lambda lambda 


1 


0. 


.4860711310 


















2 


0. 


.1073891909 


0, 


,56857 


0. 


63221 


-0. 


.23556 


-0, 


,11549 


26 


-0. 


.1229105534 


9, 


.79980 


-0. 


13175 


4. 


.66254 


4, 


.89361 


atildes 


quare <0, try 


a less extreme lambda 








? lambda 


( stop if > 


999) 


4.5 












n 


atilde2(n) 




K 


phase 


lambda lambda 


1 


0. 


.4860711310 


















2 


0. 


.1073891909 


0, 


.56857 


0. 


63221 


-0. 


.23556 


-0, 


,11549 


50000 


0. 


.2500552560 


2. 


.75176 


-0. 


28737 


4. 


.56431 


4, 


.56451 


100000 


0. 


.2499780806 


2, 


.75188 


-0. 


28663 


4. 


.56429 


4, 


.56436 


150000 


0. 


.2499892472 


2. 


.75186 


-0. 


28715 


4. 


,56432 


4, 


.56437 


200000 


0. 


.2500094861 


2. 


,75176 


-0. 


28772 


4. 


.56435 


4, 


.56439 


250000 


0. 


.2500050796 


2. 


.75188 


-0. 


28704 


4. 


.56432 


4, 


.56437 


300000 


0. 


.2499915823 


2. 


.75187 


-0. 


28722 


4. 


.56433 


4, 


.56435 


350000 


0. 


.2500079170 


2. 


.75187 


-0. 


28764 


4, 


,56435 


4, 


,56438 


400000 


0. 


.2499943758 


2, 


.75181 


-0. 


28714 


4. 


.56433 


4, 


.56435 


450000 


0. 


.2500045768 


2, 


.75184 


-0. 


28771 


4, 


,56435 


4, 


.56437 
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500000 
500001 
check: 



0.2499965741 
0.2499948277 



2.75183 -0.28722 
2.75188 -0.28741 
2.75188 -0.28769 
0.682666 B)/(A+ 



( 0.666359 A+ 

al tilde square min,max= 0.66636 
? lambda ( stop if > 999) 5 



4.56433 4.56435 

4.56434 4.56436 
4.56436 1.188E-14 3.386E-02 

1.404457 B) 



0.48607 



n 


atilde2(n) 




K 


phase 


lambda lambda 


1 


0. 


.4860711310 


















2 


0. 


.1073891909 


0. 


56857 





.63221 


-0. 


,23556 


-0. 


11549 


50000 


0. 


.2499526977 


3. 


01064 


-1, 


.08241 


4, 


.99743 


4. 


,99760 


100000 


0. 


.2500207551 


3. 


01066 


-1, 


.08330 


4. 


.99749 


4. 


,99762 


150000 


0. 


,2500144656 


3. 


01067 


-1, 


.08379 


4. 


.99752 


4. 


,99761 


200000 


0. 


,2499949927 


3. 


01057 


-1, 


.08351 


4, 


.99751 


4. 


99757 


250000 


0. 


,2499902230 


3. 


01067 


-1, 


.08422 


4. 


.99754 


4. 


99758 


300000 


0. 


,2500095958 


3. 


01063 


-1, 


.08471 


4, 


.99757 


4. 


,99760 


350000 


0. 


.2499959458 


3. 


01060 


-1, 


.08402 


4. 


.99754 


4. 


,99757 


400000 


0. 


,2500006139 


3. 


01058 


-1, 


.08470 


4. 


.99757 


4. 


,99759 


450000 


0. 


,2500018568 


3. 


01067 


-1, 


.08429 


4. 


.99755 


4. 


99758 


500000 


0. 


,2499976049 


3. 


01065 


-1, 


.08462 


4. 


.99757 


4. 


,99758 


500001 


0. 


,2500043587 


3. 


01063 


-1, 


.08476 


4. 


.99757 


4. 


,99759 


check: 






3. 


01067 


-1. 


.08480 


4. 


.99756 


7.809E-16 3.386E-02 





.666359 A+ 


0.682666 


B)/(A+ 




1.404457 B) 



( 

al tilde square min,max= 0.66636 0.48607 
? lambda ( stop if > 999) 9999 

input alpha, beta, gamma, xtilde0,nmax (stop if nmax<=0) 

STOP: 

/ul8/grpanma/magnus/jacobi@uxl2 [4] #exit 
exit 

script done on Wed Mar 2 18:01:22 1994 
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